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Abstract.  Chemical  vapor  deposition  (CVD)  process  composes  a  complex  system,  where  chemical  reaction  and  heat  and 
mass  transfer  interact  with  each  other.  And  these  macro-scale  phenomena  are  deeply  related  to  micro-scale  mechanics. 
Hence  multi-scale  analysis  is  required  to  understand  these  complicated  phenomena  and  to  develop  full-scale  simulator  of 
the  CVD  reactor.  In  this  paper,  we  present  the  macro-scale  simulation  by  the  DSMC  method.  In  those  reactors, 
sometime  the  important  species  such  as  the  reactive  intermediates  have  extremely  low  density  ratio.  This  causes  the  large 
statistical  fluctuation  in  the  DSMC  method,  where  the  number  of  particles  and  the  calculation  time  are  limited.  We 
propose  a  new  numerical  method  for  this  kind  of  problem  and  the  whole  process  of  silicon  CVD  is  simulated  by  the  new 
method.  We  simulate  the  following  CVD  process:  the  gas  mixture  of  silane  and  hydrogen  forms  a  free  expansion  jet 
through  a  nozzle  orifice  at  the  top  of  the  reactor  and  interact  with  the  heated  substrate  that  is  set  vertical  to  flow,  where 
silane  decomposes  into  silylene  and  silane  and  silylene  deposit  onto  the  surface.  It  is  confirmed  that  the  new  method  is 
very  effective  and  make  it  possible  to  analyze  the  CVD  process  more  precisely. 


1.  INTRODUCTION 

The  importance  of  controlling  and  predicting  thin  film  growth  processes  in  the  semiconductor  industry  goes  on 
increasing,  as  the  next  generation  devices  require  more  increased  component  density  and  more  decreased  feature 
size.  The  diameter  of  silicon  wafer  is  200-300  mm  and  the  thickness  of  thin  film  is  a  few  nanometers  in  the  present. 
But  it  is  said  that  the  film  thickness  becomes  less  than  1  nm  in  the  near  future.  Chemical  vapor  deposition  (CVD) 
process  is  generally  used  to  deposit  many  kinds  of  such  a  film.  Because  CVD  process  is  a  complex  system  that 
involves  chemical  reaction  and  heat  and  mass  transfer,  the  details  of  the  phenomena  in  the  reactor  are  not  understood 
well  in  spite  of  the  widespread  commercial  use.  Therefore  the  reactors  have  been  mainly  designed  empirically  until 
now.  However  it  can  be  done  more  efficiently  and  accurately  using  numerical  method  without  having  to  rely  on 
empirical  data  and  it  will  save  much  time  and  cost  to  develop  the  new  devices. 

The  studies  about  modeling  of  CVD  process  can  be  divided  broadly  into  two  groups.  One  is  the  study  of  the 
macro-scale  analysis.  This  is  about  the  whole  flow  filed  in  the  reactor,  that  is,  the  spatial  distribution  of  velocity, 
energy  and  density  from  the  viewpoint  of  reactor  size.  These  are  evaluated  by  Navier-Stokes  equations  [1H31  or  the 
DSMC  method  .  The  other  is  the  study  of  the  micro-scale  analysis  and  this  contains  intermolecular  potential 
and  collision  dynamics  from  the  viewpoint  of  molecular  size  [7,[81.  In  the  actual  CVD  process  these  physical 
phenomena  strongly  interact  each  other  and  the  multi-scale  analysis  must  be  needed  as  shown  in  Fig.l.  In  the  CVD 
reactor,  large  amount  of  molecular  collisions  result  in  the  deposition  of  thin  film.  And  the  molecular  motion  is 
controlled  by  nucleus  or  electron.  But  it  needs  enormous  computational  cost  to  solve  the  Shrodinger  equations  in 
order  to  obtain  the  whole  flow  filed  in  reactor.  Then  it  is  useful  to  extract  the  important  element  and  to  compose  the 
inter-molecular  potential.  The  same  idea  is  also  effective  for  composing  the  molecular  collision  model.  Thus  we 
can  obtain  the  physical  quantities  without  empirical  data.  Therefore  our  aim  is  to  couple  the  different  kind  of  scales 
and  develop  the  full-scale  model  to  predict  the  accurate  performance  of  the  CVD  reactor.  The  work  reported  here 
describes  macro-scale  analysis  of  low-pressure  neutral  CVD  process  using  the  DSMC  method. 

The  deposition  of  polycrystalline  silicon  from  silane  is  one  of  the  most  studied  CVD  process  and  a  large  number 
of  studies  have  been  made  by  the  DSMC  method  over  the  last  decade  or  so.  Despite  many  studies  most  of  them  are 
limited  to  the  qualitative  ones  because  the  key  gas-phase  reactions  and  the  surface  reaction  occur  infrequently  and 
the  large  differences  of  number  density  between  species  cause  statistical  problem  and  need  enormous  computational 
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cost.  In  the  process  of  polycrystalline  silicon,  the  density  ratio  of  silylene  to  silane  is  much  smaller  than  unity  in 
some  condition.  But  the  reactive  sticking  coefficient  of  silylene  is  much  larger  than  that  of  silane  and  the  difference 
of  density  is  offset.  As  to  the  contribution  to  the  deposition,  therefore  both  of  them  should  be  modeled  accurately. 
In  this  paper  we  introduce  a  new  DSMC  schemes  for  the  trace  species  chemistry  using  particle  dependent  weight 
factors  in  order  to  overcome  this  difficulty.  First,  the  new  method  is  described  and  the  differences  between  the  new 
method  and  the  conventional  one  are  discussed.  Next,  low-pressure  silicon  CVD  process  is  modeled  by  this  new 
method  and  the  effect  is  evaluated. 


FIGURE  1.  Schematic  diagram  of  the  multi-scale  analysis 


2  GAS-PHASE  REACTION  WITH  WEIGHT  FACTOR 

First,  let  us  suppose  a  typical  irreversible  dissociation  reaction  as 

AB  +  M  ->  A  +  B  +  M.  (1) 

Here,  AB  represents  the  dissociating  molecule,  A  and  B  are  atoms  or  molecules,  and  M  is  the  third  body  atom  or 
molecule.  A  basic  chemistry  model  in  DSMC  method  is  the  total  collision  energy  (TCE)  model  proposed  by  Bird  P1. 
In  this  model,  the  total  energy  Ec  is  compared  with  activation  energy  Ea.  When  Ec  is  greater  than  Ea,  the  probability 
of  reaction  Pr  is  calculated.  The  probability  Pr  is  obtained  by  integrating  the  known  equilibrium  distribution  f(EJ 
for  all  the  total  collision  energy  and  equating  it  to  a  chemical  rate  constant  k(T)  at  temperature  T  as 

k{T)=  [  algPrf(Ec)dEc  =  AT”  exp{-Ea/kT),  (2) 

where  a,  is  the  collision  cross  section,  g  is  the  relative  velocity,  Ea  is  the  activation  energy  of  the  reaction,  A  and  77 
are  constants,  and  k  is  the  Boltzmann  constant.  As  mentioned  above,  when  the  reaction  is  infrequent,  that  is, 
probability  Pr  is  low,  it  causes  the  large  difference  of  density  between  the  chemical  species.  Hence  the  species 
dependent  weight  factor  Wt  (;=AB,  A,  B,  M)  is  usually  applied  to  evade  the  statistical  problem.  In  the  following 
discussion,  we  suppose  that  species  AB  is  diluted  gas  and  species  M  is  carrier  gas  and  that  WM>  WAB>  WA=  WB.  Let  us 
start  with  the  calculation  of  the  probability  Pr  in  the  following  section. 

2.1.  Conventional  Method 


In  the  TCE  model  probability  Pr  is  assumed  to  have  the  form 


P=Q 


\n+i,j+y2 


1-- 


(3) 


where  C 1  is  constant,  which  is  determined  from  Eq.  (1)  and  Qj  is  the  average  number  of  internal  degree  of  freedom, 
which  is  arbitrary  parameter  in  this  model  and  is  the  viscosity-temperature  exponent.  In  the  conventional  method, 
probability  Pr  is  determined  by  Eq.  (3)  every  time  the  collision  with  Ec>Ea  occur  and  compared  with  the  random 
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number  to  decide  if  the  collision  is  reactive  or  not.  If  the  collision  is  reactive,  the  reactant  AB  vanishes  from  the 
computational  region  and  the  products  A  and  B  are  created.  However,  this  method  has  severe  problem  when  the 
ratio  of  weight  factor  of  species  A  and  B  to  species  AB  is  far  from  unity.  For  example,  when  WAB=1.0xl04  and 
WA=WB=1.0,  the  10,000  computational  particles  of  species  A  and  B  are  created  in  the  same  position  and  at  the  same 
time.  As  a  result,  tens  of  thousands  of  computational  particles  or  enormous  computational  time  must  be  needed  to 
collect  the  statistic.  The  engineering  DSMC  simulations  cannot  afford  to  accept  such  a  requirement. 

2.2.  New  Method 

We  propose  the  new  method  to  deal  with  the  trace  species  chemistry  much  more  efficiently.  The  basic  idea  of 
the  new  method  is  gradual  decrease  of  the  reactants  weight  every  reactive  collision,  and  this  method  consists  of  three 
steps.  In  the  first  step,  the  expected  value  that  one  computational  particle  is  involved  in  the  reaction  is  evaluated  and 
the  computational  particle  is  divided  into  two  parts.  In  the  next  step,  chemical  reaction  occurs  in  one  part  and 
another  part  is  regarded  as  elastic  or  inelastic  collision  with  the  conventional  method  by  Bird  [9).  The  last  step  is  to 
merge  the  two  parts  if  there  is  some  remainder  of  reactant  after  reactive  collision. 

2.2.1.  Stepl 

Again  we  consider  the  chemical  reaction  expressed  in  Eq.  (1).  The  first  step  is  to  divide  the  computational 
particle  of  AB  and  M  into  two  parts.  If  we  suppose  that  a  computational  particle  AB  is  a  set  of  Wr  virtual  particles, 
where  Wr=WAB/WA,  the  number  of  virtual  particles  involved  in  chemical  reaction  is  expected  from  Eq.  (3)  as, 

/  y +(«+y2 

(K)  =  WJW,  Pr  =  WrC,  1  --f  £/"•-' .  (4) 

l  Ec) 

This  means  that  <Wp>  virtual  particles  out  of  Wx  are  involved  in  reaction.  From  this  discussion  the  computational 
particle  of  AB  and  M  are  divided  into  the  virtual  two  particles  as  shown  in  the  left  side  of  Fig.  2.  One  is  the  part  that 
is  involved  in  chemical  reaction  with  the  weight  factor  expressed  as  Eq.  (5).  And  the  collision  of  another  part  is 
regarded  as  elastic  or  inelastic  and  the  weight  factor  is  expressed  as  Eq.  (6). 

K,=K=(K) 

w:b  =  wab-^,  w-=wu-(wi .} 

Though  the  computational  particles  are  divided,  the  mass  conservation  is  still  satisfied. 

2.2.2.  Step  2 

In  this  step,  collisions  are  performed  using  the  conventional  DSMC  method  for  both  parts  as  shown  in  Fig.2. 
First,  we  consider  the  reactive  collision  part.  As  mentioned  above,  virtual  particles  ABr  and  Mr  have  the  same 
weight.  Hence  the  conventional  TCE  model  can  be  applied,  and  species  A  and  B  are  created  with  the  weight  factor 
WA  and  WB.  The  number  of  the  created  computational  particles  of  A  and  B  is  1EAB.  Though  the  virtual  particle  ABr 
vanishes,  the  virtual  particle  Mr  still  remains  with  the  weight  of  JTV  Next,  nonreactive  collisions  for  another  part 
are  considered.  Because  the  weights  of  the  virtual  particles  ABn  and  Mn  are  not  the  same,  the  linear  momentum  and 
energy  conservation  are  not  satisfied  with  the  conservative  method.  The  conservative  weighting  scheme  by  Boyd  [1 11 
gives  good  solution  to  this  problem.  This  scheme  conserves  both  linear  momentum  and  energy  during  collisions. 
Therefore  all  the  post-collision  and  post-reaction  velocity  and  energy  are  determined.  Though  all  the  products  A  and 
B  are  created  without  consuming  the  reactants  AB,  the  mass  conservation  is  hold  by  decreasing  the  weight  of  AB 

2.2.3.  Step  3 

The  last  step  is  to  determine  the  velocity  and  energy  of  third  body  M  as  shown  in  Fig.2.  Two  virtual  particles  Mr 
and  Mn  have  the  velocity  vector  CM  and  C°M,  and  have  the  different  weight  WfM,  IE'm  respectively.  Again  the 


(5) 

(6) 


208 


conservative  weighting  scheme  by  Boyd  can  be  applied  to  merge  together  the  two  virtual  particles.  The  final 
velocity  vector  Cu  is  evaluated  as 


CM=(\-</>)CnM+</>CrM  (7) 

where  0=<Wp>/WM.  And  the  energy  loss  between  the  pre-  and  post-collision  state  is  written  as 

4E=£»*ut(i-t)c?cr.  (8) 

i 

This  energy  loss  is  added  by  the  same  way  as  Boyd’s  method.  And  the  internal  energy  is  also  determined  similarly. 
This  step  can  be  omitted  if  the  reaction  does  not  involve  the  third  body  like  the  exchange  reaction. 


It  is  noted  that  each  computational  particle  of  products  A  and  B  will  be  created  each  time  step  with  the  appropriate 
weights  and  that  the  exact  conservation  of  momentum  and  energy  is  satisfied  at  each  collision  by  this  new  method. 
Although  the  basic  concept  of  this  method  is  the  same  as  the  integral  detailed  balancing  concept  by  Bartel  [l01,  this 
method  is  based  on  the  collision  dynamics  and  broadly  applicable. 


^MmM  reactive 


(in)  elastic 
collision 


WnABmAB 


^ABmAB 


wm">m 


WMmM  W"MmM 


Step  1 


Step  2 


Step  3 


FIGURE  2.  Schematic  diagram  of  the  new  trace  chemistry  method 


2.3.  Application  Results 


We  have  compared  the  new  method  with  the  conventional  one  by  calculating  the  reaction  rate  constant  of  Eq.  (1). 
The  test  calculation  is  performed  in  thermal  equilibrium  single  cell  at  temperature  range  from  1000K  to  1500K.  The 
initial  component  of  gas  mixture  consist  of  is  AB  (1%)  and  species  M  (99%),  and  the  weight  factor  is  WAB:  WM=l :99. 
Initially  each  number  of  computational  particles  is  10,000  respectively.  The  number  of  species  AB  gradually 
decreases  and  the  products  A  and  B  are  created  every  time  step.  The  reaction  rate  constant  is  given  by  Eq.(2),  where 
the  values  are  set  to  ^4=1.5x10  14,  rj=  0,  £'a=3.0xl0'19.  Each  calculation  has  been  continued  until  the  calculated 
reaction  rate  reaches  almost  constant.  The  left  figure  of  Fig. 3  shows  that  a  good  agreement  with  theory  is  obtained 
for  both  methods.  It  is  clear  that  the  new  method  can  reproduce  the  molecular  collision  correctly.  And  in  the  right 
of  Fig.3,  the  new  method  gives  the  results  much  faster  than  the  conventional  method.  In  the  very  low  reaction 
probability  at  1000K  the  calculation  efficiency  increases  by  a  factor  of  100. 


FIGURE  3.  Calculation  of  the  chemical  rate  constant  and  the  velocity  distribution 
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3.  SURFACE  REACTION  WITH  WIGHT  FACTOR 


Let  us  consider  the  case  that  the  molecules  adsorb  to  the  surface  with  the  reactive  probability  y.  If  y  is  much 
smaller  than  unity,  enormous  computational  time  must  be  needed  until  the  certain  numbers  of  molecules  adsorb. 
And  when  the  weight  factor  is  used  to  evade  that  problem,  the  statistical  problem  still  remains  like  the  gas-phase 
reaction,  that  is,  the  enormous  number  of  species  absorb  at  the  same  point  and  the  same  time.  But  the  new  method 
can  be  the  solution  to  this  kind  of  problem.  Consider  the  case  that  a  computational  species  A  with  weight  WA  adsorb 
to  the  surface  with  the  probability  y.  If  we  assume  that  a  computational  particle  A  is  a  set  of  WA  virtual  particles,  the 
number  of  virtual  particles  that  adsorb  to  surface  every  time  the  particle  A  reach  the  surface  is  expected  as, 

{WA)  =  WAy.  (9) 

Then  the  computational  particle  A  is  divided  into  two  part  again.  That  is,  one  part  is  adsorbed  to  the  surface  with 
weight  <WA>,  and  another  escapes  from  the  surface  with  weight  IVA-<WA>.  Thus  with  the  appropriate  weight 
virtual  particles  will  adsorb  each  time  step  and  decreasing  the  weight  of  computational  particle  A  satisfies  the  exact 
mass  conservation. 

We  have  compared  this  method  with  the  conventional  one  by  calculating  the  distribution  of  the  number  of 
molecules  that  adsorb  to  the  surface.  The  test  calculation  is  performed  in  the  single  cell  with  the  wall  at  one  side  of 
the  boundary.  When  a  particle  reaches  to  the  wall,  the  probability  y  is  compared  with  the  random  number  to  decide 
if  the  particle  adsorb  or  not.  If  it  is  not  reactive  impact,  the  diffusion  reflection  is  applied.  The  initial  number  of 
computational  particles  is  100,000  and  the  calculation  has  been  continued  until  the  0.1  %  of  molecules  adsorb  to  the 
wall.  As  shown  in  the  left  side  of  Fig. 4,  it  is  clear  that  the  new  method  gives  the  flat  distribution  while  the 
meaningless  result  is  only  obtained  by  the  conventional  method  in  this  calculation  time.  Fig.  5  shows  the  required 
calculation  time  to  obtain  the  flat  distribution,  and  the  new  method  gives  the  excellent  efficiency  even  if  the  very 
low  y,  and  the  calculation  time  is  reduced  by  a  factor  of  10,000. 


FIGURE  4.  Comparison  of  uniformity  and  calculation  time  as  to  surface  reaction 


4  APPLICATION  TO  SILICON  CVD  PROCESS 


4.1.  Flow  Condition  and  DSMC  Method 

Using  the  new  method  as  mentioned  above,  CVD  process  is  modeled  by  the  DSMC  method.  As  shown  in  Fig  5 
the  axisymmetric  vertical  CVD  reactor  with  cold  walls  with  temperature  Tw  is  considered.  The  heated  substrate  with 
temperature  Ts  is  set  to  vertical  to  flow  at  5d  below  from  the  r-axis  and  the  diameter  of  substrate  is  ds.  The  mixture 
of  reactant  gas  silane  and  carrier  gas  hydrogen  expands  through  an  orifice  of  diameter  dn  and  forms  a  free  expansion 
jet.  The  density  ratio  of  silane  to  hydrogen  at  inlet  is  1:99  and  the  upstream  and  downstream  pressure  is  1  Torr  and 
10  mTorr,  respectively.  The  three  chemical  reactions  in  the  gas-phase  and  at  substrate  are  considered  as 

SiH  4  (g)  +  H  2  (g)  — »  SiH  2  (g)  +  2H  2  (g) ,  (10) 

SiH4(g)  — »  Si(s)  +  2H2(g),  (11) 
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(12) 


SiH 2  (g)  — »  Si(s)  +  H2  (g) , 

where  (g)  and  (s)  mean  gas  and  solid  respectively.  In  the  gas  phase  silane  decomposes  into  silylene  by  the  collision 
with  hydrogen  only.  We  use  the  kinetic  constants  of  the  Arrhenius  relation  by  Coltrin  [1]  and  ^4=  1.64x1  O' l4m3/s,  rj=  0, 
£,=3.337x  10,9J  in  Eq.(2).  At  substrate  species  react  according  to  the  Eq.(ll)  and  (12)  and  in  accordance  with 
Moffat  and  Jensen  [I2],  we  assume  the  reaction  probability  y  as 

Yuan  =  1° .  7  sou  =  A  exPi~  A  /kTs)  (13) 

where  As=531x\0'2,  £i=1.298xlO'l9J.  Tlie  silicon  atom  is  adsorbed  at  the  impact  points  and  the  hydrogen  molecules 
are  desorbed  with  Maxwell  velocity  distribution  at  the  substrate  temperature.  The  computational  domain  includes 
the  half  cross  section  of  the  axisymmetric  reactor  in  (r,  z)  coordinate  system.  Specular  reflection  is  supposed  for  all 
the  computational  particles  which  across  the  symmetric  axis.  And  vacuum  condition  is  supposed  at  the  gap  between 
the  substrate  and  the  reactor  wall,  and  no  computational  particle  enters  from  the  gap.  The  computational  domain 
consists  of  about  l.OxlO4  grid  cell  and  the  total  number  of  computational  particles  is  about  5.0x10s.  The  total 
execution  time  is  on  Alpha  PC. 

The  DSMC  calculation  follows  the  basic  statement  of  the  Bird’s  method  [9]  except  the  chemical  reaction. 
Collisions  of  computational  particles  are  simulated  using  hard  sphere  model  because  the  lack  of  the  reliable  data  of 
Si-H  system  for  other  model  like  variable  hard  sphere  model.  And  the  null  collision  method  11  is  used  in  evaluation 
of  the  collision  number.  The  Borgnakke-Larsen  statistical  model  [14)  evaluates  energy  exchange  between 
translational  and  other  internal  modes.  For  the  chemical  reaction,  the  new  and  conventional  methods  are  applied.  In 
the  new  method,  the  weight  of  every  computational  particle  can  be  different  and  the  maximum  collision  number 
must  be  evaluated  considering  the  difference  of  weight.  When  the  two  computational  particles  A,  and  B,  that  are 
different  kind  of  species  collide,  the  maximum  number  of  collision  Mmax  is  estimated  as 

’  (14) 

where  <WA,i>  means  the  expected  value  of  weight  of  the  particles  A  in  the  grid  cell,  max{A,B}  means  the  maximum 
value  of  the  two,  NA  is  the  number  of  the  computational  particles  A  in  the  grid  cell,  g  is  the  relative  velocity,  cr  is  the 
collision  cross  section,  At  is  the  time  step,  and  V  is  the  volume  of  the  grid  cell.  For  the  collisions  of  the  same  kind  of 
species,  NB  can  be  replaced  to  (A^-l  ). 


FIGURE  5.  Reactor  geometry  and  computational  domain 

4.2.  Results  and  Discussion 

We  show  one  example  of  the  calculation  of  CVD  process.  In  this  calculation,  the  temperature  of  substrate  is 
1500K  and  the  inlet  pressure  is  1  Torr.  Fig.6  indicates  the  distribution  of  number  density  and  the  averaged 
translational  temperature.  A  free  expansion  jet  is  formed  and  the  high  temperature  area  is  appeared  at  heated  surface. 
It  is  supposed  that  the  dissociation  reaction  can  occur  in  this  area.  And  as  the  right  figure  indicates,  silane  is 
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accelerated  by  lighter  hydrogen  and  large  temperature  jump  at  the  substrate  is  present.  Also  in  the  area  of  high 
energy  silane,  the  dissociation  is  expected.  We  can  obtain  the  same  result  about  the  flow  properties  by  the  new  and 
conventional  method.  But,  the  gas-phase  reaction  does  not  occur  by  the  conventional  method  because  of  the  reason 
mentioned  above.  And  only  by  the  new  method  we  can  reproduce  the  dissociation  of  silane.  Fig. 7  indicates  the 
distribution  of  number  density  of  silylene.  We  can  see  the  high-density  area  at  the  sidewall  against  our  expectation. 
This  does  not  mean  that  the  reaction  occur  near  the  cold  sidewall.  As  shown  in  the  vector  plot  of  velocity  in  Fig. 7, 
the  most  likely  explanation  is  that  the  produced  silylene  is  brought  to  sidewall  by  the  convection.  And  the  diffuse 
reflection  of  the  sidewall  causes  the  high  density  of  silylene.  But  it  is  doubtful  whether  the  very  reactive  radical  of 
silylene  does  not  adsorb  the  cold  wall  at  300K.  It  is  the  future  work  to  inquire  this  in  detail.  Next,  let  us  discuss 
about  the  deposition  distribution  at  substrate.  As  shown  in  Fig.  8,  the  distribution  of  the  deposits  from  silane  is 
extremely  uniformity.  On  the  other  hand,  the  reaction  probability  of  silylene  is  much  higher  and  the  distribution 
reflects  the  number  density  in  the  gas-phase.  The  rough  distribution  of  deposits  from  silylene  comes  from  the 
statistical  fluctuation  and  it  is  not  derived  from  the  physical  phenomena.  The  optimization  of  calculation  condition 
must  be  done.  In  this  calculation  condition  the  contribution  ratio  of  silylene  to  silane  is  about  0. 1  and  it  is  very  small 
compared  with  the  atmospheric  reactor.  This  can  be  explained  by  the  small  amount  of  the  collisions  due  to  the  low 
pressure. 


FIGURE  6.  Distribution  of  the  number  density,  temperature 


FIGURE  7.  Density  distribution  of  silylene 
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FIGURE  8.  Deposition  at  heated  substrate 


5  CONCLUSION 

We  have  introduced  the  new  trace  chemistry  method  and  applied  to  the  silicon  CVD  process.  The  new  method 
was  applied  to  both  the  gas-phase  and  the  surface  reaction  and  it  is  very  effective.  By  this  method  the  enormous 
computational  cost  can  be  reduced.  In  the  analysis  of  CVD  reactor  we  succeed  to  reproduce  the  formation  of  silylene 
in  the  gas-phase.  The  optimization  of  the  calculation  condition  and  the  further  examination  of  boundary  condition 
must  be  done.  And  in  order  to  develop  the  full-scale  model,  micro-scale  analysis  including  quantum  effect  must  be 
considered.  These  are  the  perspective  for  the  near  future. 
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